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Abstract. We present a new numerical method for accurate computations of solutions to (linear) 
one dimensional Schrodinger equations with periodic potentials. This is a prominent model in solid 
state physics where we also allow for perturbations by non-periodic potentials describing external 
electric fields. Our approach is based on the classical Bloch decomposition method which allows 
to diagonalize the periodic part of the Hamiltonian operator. Hence, the dominant effects from 
dispersion and periodic lattice potential are computed together, while the non-periodic potential 
acts only as a perturbation. Because the split-step communicator error between the periodic and 
non-periodic parts is relatively small, the step size can be chosen substantially larger than for the 
traditional splitting of the dispersion and potential operators. Indeed it is shown by the given 
examples, that our method is unconditionally stable and more efficient than the traditional split- 
step pseudo spectral schemes. To this end a particular focus is on the semiclassical regime, where 
the new algorithm naturally incorporates the adiabatic splitting of slow and fast degrees of freedom. 
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1. Introduction. One of the main problems in solid state physics is to describe 
the motion of electrons within the periodic potentials generated by the ionic cores. 
This problem has been studied from a physical, as well as from a mathematical point 
of view in, e.g., [TJ El Ell [30l [34] , resulting in a profound theoretical understanding of 
the novel dynamical features. Indeed one of the most striking effect, known as Petri's 
substitution, is a modification of the dispersion relation for Schrodinger's equation, 
where the classical energy relation Ef lee (k) — ^\k\ 2 has to be replaced by the E m (k), 
to G N, the energy corresponding to the mth Bloch band [5] . The basic idea behind this 
replacement is a separation of scales which is present in this context. More precisely 
one recognizes that experimentally imposed, and thus called external, electromagnetic 
fields typically vary on much larger spatial scales than the periodic potential generated 
by the cores. Moreover this external fields can be considered weak in comparison to 
the periodic fields of the cores [2]. 

To study this problem, consider the Schrodinger equation for the electrons in a 
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semiclas steal asymptotic scaling [T21 [301 132] , J -e- in. d = 1 dimensions 



(1.1) 




ied t ip 



t=o 



— d xx ip + V T (J) V + U(x)i), 



where < £ < 1, denotes the small semiclassical parameter describing the micro- 
scopic/macroscopic scale ratio. The (dimensionless) equation (jl.ip consequently de- 
scribes the motion of the electrons on the macroscopic scales induced by the external 
potential U(x) G M. The highly oscillating lattice-potential Vr(y) € R is assumed to 
be periodic with respect to some regular lattice T. For definiteness we shall assume 
that 



i.e. r = 2ttZ. In the following we shall assume "4>\ n 6 L 2 (M.), such that the total mass 
is Mi n = || "0m 11^2 = 1, a normalization which is henceforth preserved by the evolution. 

The mathematically precise asymptotic description of ?/>(i), solution to (11.11) . as 
e —> 0, has been intensively studied in, e.g., [T71 [5TJ [3D], relying on different an- 
alytical tools. On the other hand the numerical literature on these issues is not so 
abundant [THl El E0] . Here we shall present a novel approach to the numerical treat- 
ment of (|1.1[) relying on the classical Block decomposition method, as explained in 
more detail below. The main idea is to treat in one step the purely dispersive part 
oc d xx of the Schrodinger equation together with the periodic potential Vr, since this 
combined operator allows for some sort of "diagonalization" via the Bloch transfor- 
mation. The corresponding numerics is mainly concerned with the case e <C 1 but 
we shall also show examples for a rather large e = h. Our numerical experiments 
show that the new method converges with Ax = 0(e) and At = 0(1), the latter 
being a huge advantage in comparison with a more standard time-splitting method 
used in [18l HH [20], and which usually requires At — 0(e). Moreover we find that 
the use of only a few Bloch bands is mostly enough to achieve very high accuracy, 
even in cases where U[x) is no longer smooth. We note that our method is uncondi- 
tionally stable and comprises spectral convergence for the space discretization as well 
as second order convergence in time. The only drawback of the method is that we 
first have to compute the energy bands for a given periodic potential, although this 
is needed only in a preprocessing step rather than during the time marching. On the 
other hand, this preprocessing also handles a possible lack of regularity in Vr, which 
consequently does not lead to numerical problems during the time-evolution. In any 
case the numerical cost of this preliminary step is much smaller than the costs spend 
in computing the time-evolution and this holds true for whatever method we choose. 

We remark that linear and nonlinear evolutionary PDEs with periodic coefficients 
also arise in the study of photonic crystals, laser optics, and Bose-Einstein condensates 
in optical lattices, cf. [TUJ [121 HI] and the references given therein. We expect that 
our algorithm can adapted to these kind of problems too. Also note, that in the case 
of a so-called stratified medium, see, e.g., [6], an adaptation of our code to higher 
dimensions is very likely. Finally, the use of the Bloch transformation in problems 
of homogenization has been discussed in (TJJ [TS] and numerically studied in [T3] for 
elliptic problems. Our algorithm might be useful in similar time-dependent numerical 
homogenization problems. 

The paper is organized as follows: In Section [5J we recall in detail the Bloch- 
decomposition method and we show how to numerically calculate the corresponding 



(1.2) 



V r (y + 2n) =Vv{y) Vy £ R, 
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energy bands. Then, in Section [3] we present our new algorithm, as well as the 
usual time-splitting spectral method for Schrodinger equations. In section [U we show 
several numerical experiments, and compare both methods. Different examples of U 
and Vr are considered, including the non-smooth cases. Finally we shall also study 
a WKB type semiclassical approximation in Section [5] and compare its numerical 
solution to solution of the full problem. This section is mainly included since it gives 
a more transparent description of the Bloch transformation, at least in cases where a 
semiclassical approximation is justified. 

2. The emergence of Bloch bands . First, let us introduce some notation used 
throughout this paper, respectively recall some basic definitions used when dealing 
with periodic Schrodinger operators [2J [JJ [3H 123] ■ 
With Vr obeying (|1.2[) we have: 

• The fundamental domain of our lattice T = 2irZ, is C — (0, 2ir). 

• The dual lattice T* can then be defined as the set of all wave numbers k £ R, 
for which plane waves of the form exp(ikx) have the same periodicity as the 
potential Vr- This yields r* = Z in our case. 

• The fundamental domain of the dual lattice, i.e. the (first) Brillouin zone, 
B = C* is the set of all k £ M closer to zero than to any other dual lattice 
point. In our case, that is B = (— s, |)- 

2.1. Recapitulation of Bloch's decomposition method. One of our main 
points in all what follows is that the dynamical behavior of (jl.ip is mainly governed by 
the periodic part of the Hamiltonian, in particular for e <C 1. Thus it will be important 
to study its spectral properties. To this end consider the periodic Hamiltonian (where 
for the moment we set y — x/e for simplicity) 

(2-1) H = -±d yy + V r (y), 

which we will regard here only on L 2 (C). This is possible since due to the periodicity 
of Vr which allows to then to cover all of R by simple translations. More precisely, for 
k £ B = [ — 5) i] we equip the operator H with the following quasi-periodic boundary 
conditions 

f ^(t,y + 2ir)=e 2ik ^(t iy ) Vy £ R, k £ B, 
\d v ^(t,y + 2ir)=c 2ik7r d y iP{t,y) Vy £ M, fc £ £. 

It is well known |33) that under very mild conditions on Vr, the operator H admits 
a complete set of eigenfunctions (p m (y, k),m £ N, providing, for each fixed k £ B, an 
orthonormal basis in L 2 (C). Correspondingly there exists a countable family of real- 
valued eigenvalues which can be ordered according to E\{k) < -Ea (&) — ""' — E m {k) < 
• • • , m £ N, including the respective multiplicity. The set {E m (k) | k £ B} C R is 
called the mth energy band of the operator H and the eigenfunctions ip m (-,k) is 
usually called Bloch function. (In the following the index m £ N will always denote 
the band index.) Concerning the dependence on k £ B, it has been shown [331 that 
for any m £ N there exists a closed subset A C B such that: E m (k) is analytic and 
<p m (-, k) can be chosen to be real analytic function for all k £ B\A. Moreover 

(2.3) E m —x < E m (k) < E m+1 (k) Vfc £ B\A. 

If this condition indeed holds for all k £ B then E m {k) is called an isolated Bloch band 
[32j . Moreover, it is known that 

(2.4) meas„4 = meas{fc £ B \ E n (k) = E m {k), n ^ m} = 0. 
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In this set of measure zero one encounters so called band crossings. Note that due to 
(|2.2p we can rewrite ^ m {y, k) as 

(2.5) v m (y,k) = j ky Xm (y,k) VmeN, 

for some 27r-periodic function x m (-,fc). In terms of Xm(y,k) the Bloch eigenvalue 
problem reads 



(2.6) 



H{k)xm(y,k) =E m (k)x m (y,k), 
X m (y + 2?r, fc) =Xm{y,k) Vk&B, 



where H(k) denotes the shifted Hamiltonian 

(2.7) H{k):= l -{-id v + kf + V T {y). 

Let us know introduce the so-called Bloch transform T of some function ip(t, •) £ 
L 2 (R), for any fixed t £ R, as can be found in, e.g., [30j [32] . (Some other variants of 
this transformation can also be found in the literature.) The Bloch transformation T 
is just the regular Fourier transform T on the factor £ 2 (T) followed by a multiplication 
with e~ lyk , i.e. 

(2.8) (7V0(t, k, y) := ]T ^(t, y + 2tt 7 ) e"^ 2 ^ , yeC, keB. 
It is then easy to see that 

(2.9) TUT' 1 = H{k). 

which provides a link between the eigenvalue problem (|2 . 6[) and the periodic part of 
our Schrodinger equation acting on ip(t, •). 

Most importantly though the Bloch transformation allows to decompose our orig- 
inal Hilbert space T-L = L 2 (R) into a direct sum of, so called, band spaces, i.e. 
(2-10) 

L 2 {R) = H m , H m := | ^m(t, V)= J B /(*> fc ) Mtf. fc ) dfc > /(*' ') e i2 ( B )} » 

for any fixed t 6t. This is the well known Bloch decomposition method, which implies 
that 

(2.11) V# r )eL 2 (K): i>(t,y)=^if m {t,y), ^ m eH m - 

The corresponding projection of ip(t) onto the mth band space is thereby given as 



(2.12) ip m {t, y) = (P m .V)(i, y) = I ( / Ht, QVm (C, k) d( ) Vm {y, k) dk 
and we consequently denote by 

(2.13) C m (t,k) := / ^(i,C)^ m (C,fc)dC 



the coefficients of the Bloch decomposition. For a complete description and a rigorous 
mathematical proof of this decomposition we refer to, e.g., [31], chapter XI. Here it 
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is only important to note that the Bloch transformation allows to obtain a spectral 
decomposition of our periodic Hamiltonians H, upon solving the eigenvalue problem 
(|2.6p . Roughly speaking T can be seen as some sort of Fourier transform adapted to 
the inclusion of periodic coefficients (potentials) . 

This consequently implies that, if U = 0, we can indeed Bloch transform the whole 
evolution problem (jl.l[) and decompose it into the corresponding band spaces "H m , 
i.e. we gain some sort of "diagonalization" for our evolution problem. In this case 
each ip m (t, ■) £ W m then evolves according to the newly obtained PDE 



(2.14) 




E m (-id y )ip m , i/eK, t G 

: (P m ^in)(?y). 



Here E m (—id y ) denotes the pseudo-differential operator corresponding to the (Fourier- 
) symbol E m (k), cf. [XT] [3U1 122] • The above given evolution equation comprises a 
rigorous justification of Peirl's substitution. Moreover (|2.14j) is easily solved invoking 
the standard Fourier transformation T on L 2 (R) , which yields 

(2.15) ^ m {t, y) = T- 1 (e- iB -W*/ £ (J-(P^in))(fc)) • 

Here the energy band E m (k) is understood to be periodically extended on all of R. 
To this end, note that the following relation holds 

(2.16) HTpmW, k) = e- iB ™W*/ £ C m (0, k)(T Xm )(0, k), 

as can be shown by a lengthy but straightforward calculation. 

Of course if U ^ (the non-periodic part of the potential) the time evolution 
(jl.ip in general mixes all band spaces W m , i.e. we can no longer hope to be able 
to diagonalize the whole Hamiltonian operator (which now involves also non-periodic 
coefficients). On the other hand, since U(x) = U{ey) varies only slowly on the fast 
(periodic) scale y — x/e, one might hope that even if U ^ 0, the effective Schrodinger 
type equation 



(2.17) 



ie^Cf = E m (-id y )r* + U(ey)1>%, yeR,te 
CfL n = (P m <Ain)(2/), 



holds true, at least approximately for small s <C 1. In other words, we expect the 
slowly varying external potential to be almost constant on the lattice scale and thus 
yielding only a small perturbation of the band structure determined via (|2.1j) . Indeed 
this is the case as has been rigorously proved in [12j [211 ED] , using different analytical 
approaches, (for a broader overview, see ^32 and the references given therein), where 
it is shown that 

(2-18) sup||(P m V)W -CfW|U R) < 0(e), 

tei y ' 

holds true for any finite time-interval I C K. Here ip(t) is the solution of the full 
Schrodinger equation and ^ s t ne solution of the effective model (|2.17|) . To this 

end one hast to assume that the m'th energy band is isolated from the rest of the 
spectrum though. If this is not the case, energy transfer of order 0(1) can occur at 
band crossings, the so-called Landau-Zener phenomena. 
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2.2. Numerical computation of the Bloch bands. As a preparatory step 
for our algorithm we shall first calculate Bloch's energy bands E m (k) numerically as 
follows. Analogously to []j3[27], we consider the potential Vr € C 1 (IR) and expand it 
in its Fourier series, i.e. 

(2.19) V r (y) = ]T F(A) e iA ^, V(X) = ±- [ * V r (y) dy. 

27T In 

AeZ JU 

Likewise, we expand any Bloch eigenfunctions Xm{'> k), in its respective Fourier series 



1 



2tt 



(2.20) Xm(y,k) = J2x m (\,k)e iX y, £ m (A,fc) = ^/ Xm (y , k) dy . 

(The latter should not be confused with the so-called Wannier functions which are 
given as the Fourier transformation of tp m w.r.t to k G £>.) Clearly the Fourier 
approximation of Vr, and thus also the one of Xm, depends on the regularity of Vr- If 
Vr S C°°(R) the corresponding Fourier coefficients V(X) decay faster than any power, 
as A — ► ±oo, and thus we only need to take into account a few coefficients in this case. 

For A G {—A, • • • , A — 1} C Z, we consequently aim to approximate the Sturm- 
Liouville problem (|2.6[) . by the following algebraic eigenvalue problem 



(2.21) 



H(fc) 



/ Xm(-A) \ 
Xm(l-A) 



V Xm(A-l) J 



= E m (k) 



( Xm(-A) \ 
Xm(l-A) 



V Xm(A-l) J 



where the 2A x 2A matrix H(fc) is given by 
(2.22) 

/ V>(0)+J(fc-A) 2 V(-l) 



H(fc) 



V(l) 



\ V(2A-1) 



V(0) + i(fc-A + l) 5 



V(2A-2) 



V(l - 2A) 
9(2 - 2A) 

T/^ + ^fc + A-l) 2 / 



The above given matrix H(fc) comprises 2A eigenvalues. Clearly, this number has to 
be large enough such that all the eigenvalues E m (k) which we need to use in our simu- 
lations below are counted, i.e. we need m < 2A. The numerical cost for this algebraic 
problem is about C(A 3 ), cf. [23] ■ Note however that this is the most expensive case, 
which becomes considerably smaller if one exploits possible symmetries within the po- 
tential Vr, cf. Example 14.11 below (see also [TU1 [23 1221 ISS] ) • In any case the number A 
is independent of the spatial grid, thus the numerical costs of this eigenvalue problem 
are almost negligible compared to those spend in the evolutionary algorithms below. 
The approximate numerical computations of the Bloch bands E m (k) can be seen as 
a preprocessing, to be done only once and remain unchanged as time evolves. 

Remark 2.1. Accurate computations of the energy bands needed in practical 
applications, i.e. in more than one spatial dimensions and for different kind of 
(composite) material, becomes a highly nontrivial task. Nowadays though, there al- 
ready exists a huge amount of numerical data comprising the energy band structure 
of the most important materials used in, e.g., the design of semiconductor devices, 
cf. fIR Wb\ We note that some of these data is available online via the URL 



http : //www. research. ibm. com/DAMOCLES/home .html or http : / / cmt . dur . ac .uk/sjc 
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d also http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html. In the context of 



an 

photonic crystals the situation is similar J22f . Thus, relying on such data one can 
in principle avoid the above given eigenvalue- computations (and its generalizations 
to more dimensions) completely. To this end, one should also note that, given the 
energy bands E m (k), we do not need any knowledge about Vy in order to solve (|l.ll) 
numerically, cf. the algorithm described below. 

3. Bloch decomposition based algorithm vs. time-splitting spectral 
methods. For the convenience of computations, we shall consider the equation II) 
on a bounded domain T>, say on the interval T> = [— «i,«2], f° r some large enough 
K\,K2 > 0. Moreover we shall equip T> with periodic boundary conditions. However, 
this periodic computational domain T> should not be confused with the periodic struc- 
ture induced by the lattice potential. Without loss of any generality, we assume that 
V = [0,2tt]. 

For practical reasons we shall now introduce, for any fixed t 6 R, a new unitary 
transformation of ip(t, •) G L 2 (M) 

(3.1) #(t, y, k) := £ (y + 2^7)) e^ 7 , y G C, k G B, 

7GZ 

which has the properties that ip is quasi-periodic w.r.t y G T and periodic w.r.t. 
k G T*, i.e. 

(3.2) ^(t,y + 27r,fc) = e i2 ^^,y,fc), $(t,y,k + l) = $(t,y,k). 

One should note that ip is not the standard Bloch transformation T, as defined in 
(|2.8|) . but it is indeed closely related to it via 

(3.3) {Tm,y^) = i>{t,y,k)e- iyk , keB, 
for e = 1. Furthermore, we have the following inversion formula 

(3.4) ^(i,e(y + 27r 7 ))= / ^(t,y,k)e i2 ^dk, 

Jb 

which is again very similar to the one of the standard Bloch transformation [32] ■ The 
main advantage in using tp, instead of Tip itself, is that we can rely on a standard fast 
Fourier transform (FFT) in the numerical algorithm below. If one aims to use Tip di- 
rectly one would be forced to modify a given FFT code accordingly. A straightforward 
computation then shows that 

(3.5) C m {t, k) = f £(t, C, k)Tp m (C, k) dC, 

Jc 

where C m (t,k) is the Bloch coefficient, defined in (|2.13l) . 

In what follows, let the time step be At = T/N, for some N G N, T > 0. Suppose 
that there are L G N lattice cells within the computational domain V = [0, 2ir]. In 
this domain, the wave function ip is numerically computed at L x R grid points, for 
some R £ N. In other words we assume that there are R grid points in each lattice 
cell, which yields the following discretization 



(3.6) 



1 £ - 1 

ki = - - H — , where I = {1, • • • , L} C N, 

Z ij 

2it(r — 1) 

Vr = — -5 — where r = {l,--- ,R} C N, 
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and thus we finally we evaluate ip n = ip{t n ) at the grid points x = e(2irj + y), i.e. 
(3.7) xt, r = e(2ir(l-l) + y r ). 

We remark that in our numerical computations we can use R <C L, whenever e<1, 
i.e. we only use a few grid points within each cell. Now we shall describe precisely 
the Bloch decomposition based algorithm used to solve (II. ip . 

3.1. The Bloch decomposition based algorithm (BD). Suppose that at 
the time t n we are given ip(t n ,X£, r ) « 4>i r - Then ipg^ , i.e. the solution at the (next) 
time step t n+ i = t n + At, is obtained as follows: 
Step 1. First, we solve the equation 

(3-8) ieft^-y^ + Vr^)^, 

on a fixed time-interval At. To this end we shall heavily use the Bloch-decomposition 
method, see below. 

Step 2. In a second step, solve the ordinary differential equation (ODE) 

(3.9) iedtip = U{x)ij}, 

on the same time-interval, where the solution obtained in Step 1 serves as initial 
condition for Step 2. We easily obtain the exact solution for this linear ODE by 

(3.10) i){t,x) = iP{Q,x)e- iU{x)t/e . 



Remark 3.1. Clearly, the algorithm given above is first order in time. But we 
could easily obtain also a second order scheme by the Strang splitting method, which 
means that we use Step 1 with time-step At/2, then Step 2 with time-step At, and 
finally integrate Step 1 again with At/2. Note that in both cases the scheme conserves 
the particle density p(t,x) := \ip(t,x)\ , also on the fully discrete level. Indeed Step 
1 consists of several intermediate steps which we shall present in what follows: 

Step 1.1. We first compute ip at time t n by 

L 

(3.H) i=E^ e ^'' 



Step 1.2. Next, we compute the mth band Bloch coefficient C m (t, k), at time t n , 
via (13.51). i.e. 



2ir R ~ 

C m (t n , hi) « CZ, t - E Wr&r, h) e- ifc «" 
r=l 

( 3 - 12 ) R R/2-1 

r=l \=-R/2 

where for the second line we simply inserted the Fourier expansion of Xmi given in 
(|2.20[) . Note that in total we have R Fourier coefficients for Xm- Clearly this implies 
that we need A > R/2 to hold, where A is the number of Fourier modes required in 
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the numerical approximation of Bloch's eigenvalue problem as discussed in Section 
12.21 Here we only take the R lowest frequency Fourier coefficients. 

Step 1.3. The obtained Bloch coefficients are then evolved up to the time t n+1 , 
according to the explicit solution formula (|2.15p . taking into account (I2.16|) . This 
yields 

(3.13) C^+l = C^t e - iE ^ At /c. 

Step 1.4. From here, we consequently compute ip at the new time t n+1 by 
summing up all band contributions and using the analytical formulas (I2.12[) and (I2.13[) , 
i.e. 

M M R/2-1 

(3.14) C^E^^E^!/ E X m (\k e )e^+^. 

m=l m=l X=-R/2 

Step 1.5. Finally we numerically perform the inverse transformation to (I3.1[) . i.e. 
we compute ip^ 1 from ^"^T 1 . Thus from (I3.4[) . we get 

(3.15) VC 1 = \ E^ 1 ^^' 1 - 

Note that in the BD algorithm, the main numerical costs are introduced via the 
FFT in Steps 1.1 and 1.5. This also implies that on the same spatial grid, the 
numerical costs of our Bloch transform based algorithm is of the same order as the 
classical time-splitting spectral method below. Moreover, we want to stress the fact 
that if there is no external potential, i.e. U(x) = 0, then the above given algorithm 
numerically computes the exact solution of the evolutionary problem ([l.ip , which can 
be seen analogous to a standard spectral method, adapted to periodic potentials. In 
particular this fact allows us to solve the Schrodinger equation (|1.1|) for very long time 
steps, even if e is small (see the results given below). Moreover, one should note that 
a possible lack of regularity in Vr only requires numerical care when approximating 
(|2.6p by the algebraic problem (|2.21[) . In particular, Vr itself does not enter in the 
time-evolution but only E m {k). 

3.2. A simple time-splitting spectral method (TS). Ignoring for a moment 
the additional structure provided by the periodic potential Vr , one might wish to solve 
(|1.1[) by using a classical time-splitting spectral scheme. Such schemes already proved 
to be successful in similar circumstances, see, e.g., [3j HJ [19l [24] . For the purpose of a 
detailed comparison, we present this method here: 

Step 1. In the first step we solve the equation 

(3-16) iec^ = -j d xx i>, 

on a fixed time interval At, relying on the pseudo-spectral method. 

Step 2. Then, in a second step, we solve the ordinary differential equation 

(3.17) ied t i> = (Vr (^) + U{xj) ^ 

on the same time-interval, where the solution obtained in step 1 serves as initial 
condition for step 2. Again it is easily seen, that such a scheme conserves the particle 
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density. It is clear however that, due to the inclusion of Vt (-), the exact solution of 

(3.18) r/i(t, X) = fln n (x) e -Wr(x/e)+U(x))t/^ 

involves high oscillations on different length- and time-scales as e — > (which one has 
to resolve), in contrast to ([3.10)1 . where only i/e-oscillations are present. 

Remark 3.2. In our BD algorithm, we compute the dominant effects from dis- 
persion and periodic lattice potential in one step, and treat the non-periodic potential 
as a perturbation. Because the split-step communicator error between the periodic and 
non-periodic parts is relatively small, the step size can be chosen considerably larger 
than for the SP algorithm. 

Remark 3.3. Clearly, if there is no lattice potential, i.e. Vr(y) = 0, the BD 
algorithm simplifies to the described time- splitting method TS. Moreover, a second 
order second order scheme (based on the Strang splitting algorithm) can be analogously 
obtained to the one described above, see Remark \3.1[ and a comparison of these second 
order schemes gives similar results as those shown in the following. 

Remark 3.4. For the BD algorithm, the complexities of Step 1.1 and 1.5 are 
0(RL\og(Lj), the complexities of Step 1.2 and 1.4 are 0(M LR\og(R)) , and for Step 
1.3 we have 0{ML). Also the complexity of the eigenvalue problem ([2.21jl is C(A 3 ). 
However, since A (or R) is independent of e and since we only need to solve the 
eigenvalue problem (j2.21[) once in a preparatory step, the computation costs for this 
problem are negligible. On the other hand, for the TS algorithm, the complexities of 
Step 1 and 2 are 0{RL\og{RL)) and O(RL) respectively. As M and R are indepen- 
dent of e, we can use R <C L and M -C L, whenever s <C 1. Finally the complexities 
of the BD and TS algorithm in each time step are comparable. 

4. Numerical experiments. In this section, we shall use several numerical 
examples to show the efficiency of our algorithm. We shall choose for (jl.ip initial 
data V'in € S (R) of the following form 

(4.i) Mx) = (^y f \-^-«r t 

Let us perform a decomposition of -0in in terms of the Bloch bands, and take a 
summation of the first m — 1, • • • , Mq energy bands, for some finite (cut-off) number 
M £ N. A picture of the corresponding band densities p e m := P^^inl 2 is given in 
Figure SHI for m = 1, • • • ,4. Here (F^ n if>i n )(x) is the e-scaled projection onto rl e m , 
obtained from ([2.12jl by replacing ip m (y,k) — > e~ l ' 2 (p m {x/£,k). Since ipi n is smooth 
we expect that only very few bands have to be taken into account in the Bloch 
decomposition. Indeed we observe that the amount of mass corresponding to P^^n, 
i.e. the mass concentration in each Bloch band, decays rapidly as m -> oo, see Table 
14.11 In other words, the number M is essentially determined by the regularity of i/'in 
in each cell. Note that M is independent of e. 

To compute the evolution of these initial data we shall take into account M > Mq 
bands. Note that only in cases where U(x) = one can take M to be identical to 
Mq, the initial band cut-off. The reason is that if U(x) is nonzero Step 2 in the 
BD algorithm given above mixes all bands. In particular all the ip m (t) are no longer 
orthogonal to each other. Roughly speaking however, if e is very small, all band 
spaces TL m remain "almost orthogonal" and thus the mass within each Bloch band, 

II 1 1 2 

i.e. M^(<) := || P^/)(t)|| i . 2 „ t , is "almost conserved" . More precisely it is conserved up 
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111= 1 in — .2 mi— * in I 



Fig. 4.1. |P^i/> in | 2 , m =!,-•• ,4 fore = i 



Table 4.1 

The traiues of := || P^Vin Hi^py , /ore = ±: 



m 


1 


2 


3 


4 




M £ 


7.91E - 1 


1.11E - 1 


5.92E - 1 


8.80E 


- 2 


?Tt 


5 


6 


7 


8 






8.67E - 2 


2.81E - 3 


2.80E - 3 


4.98E 


- 5 



to errors 0(e) on time scales 0(1). Thus, by checking mass conservation after each 
time step one gets a rather reliable measure on the amount of mixing of the bands. 
In other words if the mass conservation after some time steps gets worse, one has to 
take into account more bands to proceed. 

We find numerically that the use of M = Mo ~ 8 bands already yields satisfactory 
results for e = In the following though we shall even compute M — 32 energy 
bands, which is by far sufficient for our purposes (even if e = |). Note that the 
number of required bands M depends on the regularity properties of U(x), as well 
as on the considered time-scales (which might be even longer than 0(1), the case 
considered here). This approximation problem is more or less analogous to the one 
appearing in spectral schemes for PDEs with non-smooth coefficients. 

Concerning slowly varying, external potentials U, we shall choose, on the one 
hand, smooth functions which are either of the form 

(4.2) U(x) = £x, 

modelling a constant (electric) force field £ € K, or given by a harmonic oscillator 
type potential 

(4.3) U(x) = \x - tt| 2 . 

On the other hand, we shall also consider the case of an external (non-smooth) step 
potential, i.e. 

(4.4) "(*) = {;; ds G e. [f ' f] 

Within the setting described above, we shall focus on two particular choices for the 
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lattice potential, namely: 

Example 4.1 (Mathieu's model). The so-called Mathieu's model, i.e. 

(4.5) V r (x) = cos(x), 

as already considered in \19^ . (For applications in solid state physics this is rather 
unrealistic, however it fits quite good with experiments on Bose-Einstein condensates 
in optical lattices.) In this case all Fourier coefficients V(\), appearing in (|2.19l) are 
zero, except for V (±1) = | and thus H(fc), given in (|2.22[) , simplifies to a tri-diagonal 
matrix. 

Example 4.2 (Kronig-Penney's model) . The so-called Kronig-Penney's model, 

i.e. 

( 4 - 6 ) V r (x) = l-J2Ke[f + 2^+2^ 

where 1q denotes the characteristic function of a set Q C R. In contrast to Mathieu's 
model this case comprises a non-smooth lattice potential. The corresponding Bloch 
eigenvalue problem is known to be explicitly solvable (see, e.g., \19f ). 

In order to compare the different numerical algorithms we denote by ip ts (t,x) the 
solution gained from the time-splitting spectral method, whereas ip hd (t, x) denotes the 
solution obtained via the new method base on Bloch's decomposition. Both methods 
will be compared to the "exact" solution ip°*(t, x), which is obtained using a very fine 
spatial grid. We consider the following errors 



\L° 



A b J/ ts (t) :=\\r x (t,-)-^ hd/ts (t,-)\ 
A h 2 d ' ts (t):=\\r x (t,-)-^ d/ts (t r) \\ L2i 



between the "exact solution" and the corresponding solutions obtained via the Bloch 
decomposition based algorithm resp. the classical time splitting spectral method. The 
numerical experiments are now done in a series of three different settings: 

• First we shall study both cases of Vr, imposing additionally U(x) = 0, i.e. no 
external potential. The obtained results are given in Table |4~51 where e = |, 

and yjj24, respectively. In the last case the oscillations are extremely spu- 
rious. As discussed before, we can use only one step in time to obtain the 
numerical solution, because the Bloch-decomposition method indeed is "ex- 
act" in this case (independently of e). Thus, even if we would refine the time 
steps in the BD algorithm we would not get more accurate approximations. 
On the other hand, by using the usual time-splitting method, one has to refine 
the time steps (depending on e) as well as the mesh size in order to achieve 
the same accuracy. More precisely we find that At = 0(e), Ax = 0(e a ), for 
some a > 1, is needed when using TS (see also the computations given in 
|19jh In particular a > 1 is required for the case of a non-smooth lattice po- 
tential Vr- (Note that if Vr = it is well known that Ax = 0(e), is sufficient, 
cf. Efl H |23|). 

• In a second series of numerical experiments we shall consider only Example 
14.11 for the periodic potential but taking into account all three cases of the 
external potentials U, as given above. In Fig. I6.1H6.6[ we show the obtained 
numerical results for e = ^, and e — j^a, respectively. We observe that, 
if e = 0(1), the Bloch-decomposition method gives almost the same results 
as time-splitting spectral method. However, if e <C 1, we can achieve quite 
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good accuracy by using the Bloch-decomposition method with At — 0(1) 
and Ax = O(e). On the other hand, using the standard TS algorithm, we 
again have to rely on much finer spatial grids and time steps to achieve the 
same accuracy. 

• We finally show the numerical results obtained by combining external fields 
and a non-smooth lattice potential given by Example 14.21 As before we 
include all three cases for the external potential U. The cases e = \ , and 
are studied and the obtained results are given in Fig. 16.7116.121 respectively. 
We observe that the results of the Bloch-decomposition are much better than 
the time-splitting spectral method, even if e = \. Moreover, as e gets smaller, 
the advantages of the Bloch-decomposition method are even better visible. 

To convince ourselves that only a few Bloch bands contribute to ||V'I|l 2 (r)j even after 
time steps 0(1), we show in the following table the numerical values of yi e m (t) = 

ii i i 2 

|| ^ > mV'(*)|| i 2( R )j f° r ni = 1,' • • ,8, corresponding to the solution of Example 14.11 with 
U given by (|jT3"l) . 

Table 4.2 

The mass ofip(t,x), solution to Examp le \4- l\ with external potential 1 14.311 . decomposed into the 
Bloch bands for e = J?; at time t = 1: 



m 


1 


2 


3 


4 




7.89E - 1 


1.10E - 2 


5.92E - 1 


9.38E - 2 


m 


5 


6 


7 


8 


M e 

m 


7.15E - 2 


3.50E - 3 


1.80E - 3 


5.63E - 5 



We also check the conservation of the total (discrete) mass, i.e. IVKOIIz^d)- We 
find that numerically it is of the order 10 -6 for the smooth lattice potential (|4.5j) and 
1CT 3 for the non-smooth case (|4.6[) . The latter however can be improved by using a 
refined spatial grid and more time steps. 

In summary we find (at least for our one dimensional computations) that, relying 
on the new Bloch-decomposition based algorithm, one can use much larger time steps, 
and sometimes even a coarser spatial grid, to achieve the same accuracy as for the 
usual time-splitting spectral method. This is particularly visible in cases, where the 
lattice potential is non longer smooth and e -C 1. Indeed in these cases the BD 
algorithm turns out to be considerably faster than the TS method. 

Remark 4.1. In view of our results the earlier numerical studies based on TS 
methods U9\ [iffi \2U$, should be taken with some care, in particular when comparing 
the full Schrodinger solution to the semiclassical approximation beyond caustics. 

5. Asymptotic analysis in the semiclassical regime. For completeness we 
shall also compare the numerical solution of the Schrodinger equation (jl.lj) with its 
semiclassical asymptotic description. To this end we shall rely on a multiple scales 
WKB-type expansion methods, even though there are currently more advanced tools 
at hand, cf. [T71 1301 122 ■ The WKB method however has the advantage of given a 
rather simple and transparent description of ip(t), solution to (|1.1[) . for e -C 1, (at 
least locally in-time). Since the Bloch decomposition method itself is rather abstract 
we include this approximative description here too, so that the reader gets a better 
feeling for the appearing quantities. Moreover this two-scale WKB method can also be 
used for nonlinear Schrodinger dynamics [12], a problem we shall study numerically 
in an upcoming work. 
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Table 4.3 

The results of Example \4 H with U(x) = 0.' 



Spatial discretization error test at time t = 1.0 for e = 1/2. 
For TS At = 0.0001 and for BD At = 1. 



mesh size Ax/e 


1/2 1/4 1/8 1/16 




^,At(*.0-^ eX (*--) 


i 2 


4.33E-1 2.53E-1 2.80E-2 6.42E-6 


convergence order 


0.8 3.2 12.1 




V-aIa^-)-^,-) 


i 2 


3.01E-1 1.95E-1 1.39E-2 1.17E-6 


convergence order 


0.6 3.8 13.5 


Spatial discretization error test at time t = 0.1 for e = 1/32. 
For TS At = 0.00001 and for BD At = 0.1. 


mesh size Ax/e 


1/2 1/4 1/8 1/16 




^,A t (*.-)-^ cx (i,') 


i 2 


2.88E-1 1.08E-1 9.63E-4 1.33E-7 


convergence order 


1.4 6.8 12.8 




v&U*(*. -)-v> ex (v) 


i 2 


2.53E-1 7.34E-2 8.97E-4 4.95E-10 


convergence order 


1.8 6.4 20.8 


Spatial discretization error test at time t = 0.01 for e = 1/1024. 
For TS At = 0.000001 and for BD At = 0.01. 


mesh size Ax/e 


1/2 1/4 1/8 1/16 




V>^,A t (t.-)-V> cx (t,-) 


i 2 


5.14E-1 1.94E-1 1.08E-3 6.08E-8 


convergence order 


1.4 7.5 14.1 




^AlAi(V)-V>° X (t,-) 


i 2 


2.64E-1 6.83E-2 2.29E-4 1.71E-10 


convergence order 


2.0 8.2 20.4 



5.1. The WKB formalism. To this end let us suppose that the initial condition 
is ol (two-scale) WKB-type. More precisely assume 

M 

(5.1) V>in(aO = $> m (*,~) e^)/ £ , 

m— 1 

with some given real-valued phase <p S C°°(R) and some given initial (complex- valued) 
band- amplitudes u m {x,y + 2ir) = u m (x,y), each of which admits an asymptotic de- 
scription of the following form 

(5.2) u m {x,y) ~u m (x,y)+eul n (x,y) + O(e 2 ) VmeN. 

Here and in the following we shall only be concerned with the leading order asymptotic 
description. 

Remark 5.1. Note that we do consider only a single initial WKB-phase <j)(x) 
for all bands m G N. We could of course also allow for more general cases, like 
one WKB-phase for each band or even a superposition of WKB-states within each 
band. However in order to keep the presentation clean we hesitate to do so. The 
standard WKB approximation, for non-periodic problems, involves real-valued ampli- 
tudes u°(x),u 1 (x), . . . which only depend on the slow scale. It is well known then, 
cf. [IHHT], that the leading order term u^, m £ N, can be decomposed as 

(5-3) u° m (x, -) = f m (x)xm (p^O)) , 
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Table 4.4 

The results of Example \4 H with linear external potential 14.21 : 



Spatial discretization error test at time t = 0.1 for e = 1/2. 
For TS At = 0.0001 and for BD At = 0.01. 



mesh size Ax/e 


1/2 1/4 1/8 1/16 




^g«,At(*.-)-r"(t,0 


i 2 


2.73E-1 9.22E-2 5.78E-3 4.73E-6 


convergence order 


1.6 4.0 10.3 




^Ai,At(*.-)-^ ex (*r) 


! 2 


3.15E-1 1.55E-1 1.32E-2 3.36E-6 


convergence order 


1.0 3.6 11.9 



Spatial discretization error test at time t = 0.01 for e = 1/1024. 
For TS At = 0.00001 and for BD At = 0.001. 



mesh size Ax/e 


1/2 


1/4 


1/8 


1/16 








5.22E-1 


1.98E-1 


1.53E-2 


3.19E-5 


convergence order 




1.4 


3.7 


8.9 




^i,At(V)-iP*(V) 


! 2 


4.71E-1 


1.61E-1 


9.17E-3 


6.08E-6 


convergence order 




1.5 


4.1 


10.6 



Temporal discretization error test at t = 0.1 for e = 1/2 and Ax/e = 1/128. 



time step At 


1/10 


1/20 


1/40 


1/80 








2.59E-4 


6.47E-5 


1.62E-5 


4.04E-6 


convergence order 




2.0 


2.0 


2.0 






! 2 


4.86E-5 


1.23E-5 


3.08E-6 


7.60E-7 


convergence order 




2.0 


2.0 


2.0 



Temporal discretization error test at t = 0.01 for e = 1/1024 and Ax/e = 1/128. 



time step At 


1/1000 1/2000 1/4000 1/8000 






! 2 


6.60E-2 1.54E-2 3.81E-3 9.45E-4 


convergence order 


2.1 2.0 2.0 


time step At 


1/100 1/200 1/400 1/800 




^AlAt(*.')-^ ex (*r) 


; 2 


3.32E-3 7.54E-4 1.42E-4 3.16E-5 


convergence order 


2.1 2.4 2.2 



where we assume the m-th energy band to be non- degenerated (for simplicity) and 
isolated from the rest of spectrum. We can choose an arbitrary f m 6 <S(R). In other 
words, there is an adiabatic decoupling between the slow scale x and fast scale x/e. 
Indeed, a lengthy calculation, invoking the classical stationary phase argument, cf. 
chapter 4.7 in !7|, shows that in this case the band projection P^ip can be approxi- 
mated via 

(5.4) P^Mx) - f m (x) Xm ^ {x)/e + 0(e). 

This approximate formula shows the origin of the high oscillations induced either 
by Vr, described by Xmi or by the dispersion, described by 4>(x). We note that in 
general the higher order terms (in e), such as u x m etc., are of a more complicated 
structure than (|5.3[) . but we shall neglect these terms in what follows (see, e.g., [12] 
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Table 4.5 

The results of Example \4-S\ with harmonic external potential 114.311 ; 



Spatial discretization error test at time t = 0.1 for e = 1/2. 
For TS At = 0.0001 and for BD At = 0.01. 



mesh size Ax/e 


1/2 1/4 1/8 1/16 




^g«,At(*.-)-r"(t,0 


i 2 


2.71E-1 8.87E-2 5.19E-3 1.32E-4 


convergence order 


1.6 4.1 5.3 




^Ai,At(*.-)-^ ex (*r) 


! 2 


3.23E-1 9.08E-2 7.03E-3 1.27E-4 


convergence order 


1.8 3.7 5.8 



Spatial discretization error test at time t = 0.01 for e = 1/1024. 
For TS At = 0.00001 and for BD At = 0.001. 



mesh size Ax/e 


1/2 


1/4 


1/8 


1/16 








3.99E-1 


3.67E-1 


2.19E-1 


1.10E-1 


convergence order 




0.1 


0.7 


1.0 




^i,At(V)-iP*(V) 


! 2 


2.06E-1 


5.64E-2 


8.16E-3 


6.40E-4 


convergence order 




1.9 


2.8 


3.7 



Temporal discretization error test at t = 0.1 for e = 1/2 and Ax/e = 1/128. 



time step At 


1/10 1/20 1/40 1/80 








1.02E-3 6.41E-4 3.80E-4 2.18E-4 


convergence order 


0.7 0.8 0.8 






! 2 


4.20E-6 1.02E-6 2.22E-7 5.56E-8 


convergence order 


2.0 2.2 2.0 



Temporal discretization error test at t = 0.01 for e = 1/1024 and Ax/e = 1/128. 



time step At 


1/1000 


1/2000 


1/4000 


1/8000 






! 2 


1.21E-1 


1.18E-1 


1.10E-1 


1.10E-1 


convergence order 




0.04 


0.1 


0.0 


time step At 


1/100 


1/200 


1/400 


1/800 




^AlAt(*.')-^ ex (*r) 


; 2 


3.30E-5 


5.21E-6 


1.23E-6 


3.16E-7 


convergence order 




2.6 


2.1 


2.0 



for more details). One consequently finds that ip(t) obeys a leading order asymptotic 
description of the form 

M 

(5.5) VM~ ^« m (i,x)Xm(^ 1 i(M))e'*"' M/E + 0(£), 

m— 1 

where 4> m {t,x) 6 C°°([0,t c ) x R) satisfies the mth band Hamilton- J acobi equation 




dt<f> m {t, x) + E m (d x (j) m ) + U{x) =0, m E N, 
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Also, the (complex-valued) leading order WKB-amplitude a m (t,x) € C°°([0,z; c ) x 1 

satisfies the following semiclassical transport equations 

(5.7) 

<9 t a m + dkE m (d x 4>m)d x a m + - d x (dkE m (d x (j) m ))a m - (f3 m (t, x)d x U(x)) a m = 0, 



I t=o 



fm(x). 



with j3 m (t,x) := (Xrn(y,k),d k Xrn(y,k)) L2 (Q, evaluated at k = d x (/) m , the so-called 
Berry phase term. 

Remark 5.2. Note that the Berry term is purely imaginary, i.e. f3 m (t,x) £ iR, 
which implies the following conservation law 

(5.8) d t \a m \ 2 + d x (d k E m {d x m )\a m \ 2 ) =0 VmsN. 

Of course the above given WKB-type expansion method is only valid up to the (in 
general finite) time < t c < oo, the caustic onset-time in the solution of (15.61) . Here 
we shall simply assume that t c > holds, i.e. no caustic is formed at time t — 0, which 
is very well possible in general. We note that in the considered numerical examples 
below we indeed have t c > and we refer to pj] for a broader discussion on this. For 
t > t c one would need to superimpose several WKB-type solutions corresponding to 
the multi- valued solutions of the flow map (x,t) M> X t (x) = X t (x;d x (j)(x)), where 



(5.9) 



X t = d k E m (^ f ), X Q = x, 

S t = - d x U(X t ), S = d x (f>(x) 



Numerically we shall use the relaxation method introduced in [3S] to solve the Hamilton- 
Jacobi equation (|5.6p . Consequently we can solve the system of transport equations 
(I5.7|) by a time-splitting spectral scheme similar to the ones used above. 



5.2. Numerical examples. We shall finally study the WKB approach, briefly 
described above, by some numerical examples. Denote by 

M 

(5.10) r c (t,x) := ]T fm{t,x) X m (p^0 m ) e *»(*.»>/«, 

m— 1 

the approximate semiclassical solution to the Schrddinger equation (jl.ip . In the fol- 
lowing examples we only take into account a harmonic external potential of the form 

(S3!). 

Example 5.1 (Mathieu's model). We first consider Mathieu's model (14.5[) and 
choose initial condition in the form 

(5.11) Vi„(aO=e- 5 (*-^ 2 Xl (^,o), 

i.e. we choose 4>{x) = and restrict ourselves to the case of only one band with index 
m = 1. (Since E\(k) is an isolated band the analytical results of ^ I VA \21f . then 
imply that we can neglect the contributions from all other bands m > 1 up to errors 
of order 0(e) in L 2 (M.) n L°°(R), uniformly on compact time-intervals.) In this case, 
we numerically find that no caustic is formed within the solution of (|5.6j) at least up 
to t = 1, the largest time in our computation. Note that (|5.11[) concentrates at the 
minimum of the first Bloch band, where it is known that 

(5.12) E m (k) « ^L + E m (0), 

2m 
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Fig. 5.1. The graph of d x <j> 2 (t, x) at t = 0.24. 

This is the so-called parabolic band approximation, yielding an effective mass m* S K. 
In Table \5.1l we show the results with an additional harmonic external potential, cf. 
(|4.3|) , for e = ^ and e = respectively. 

Table 5.1 

Difference between the asymptotic solution and the Schrddinger equation for examvle \5.1\ (At = 
1CT 4 , Ax = l/32768j: 







e 


l 

32 


l 

1024 


sup 

0<t<l 


1 il>(t, x) 


-j/j sc {t,x) \\ L2(W) 


6.68E - 3 


3.08E - 4 


sup 

0<t<l 


ip(t, x) 


-iP BC (t,x) \\ L oo (M) 


5.57E - 2 


2.38E - 3 



Note that these numerical experiments, together with those given below, confirm the 
analytical results given in [7| \1S^ \21f . 

Example 5.2 (Kronig- Penney 's model). Here, we consider again the Kronig- 
Penney's model (|4.6|) . First we use the same initial condition as given in (|5.11[) but 
with m — 2, which again corresponds to an isolated energy band. The corresponding 
numerical results for e = ^ and are s ^ own ^ n Table 

Table 5.2 

Difference between the asymptotic solution and the Schrddinger equation for examvle \5.2\ for 
initial condition l5~TH (At = 1CT 4 , Ax = 1/32768J: 





£ 




1 

32 


l 

1024 


sup 

0<i<0.1 


1 il>(t, x) - 


4> sc (t,x) \\ L2{R) 


1.18E - 2 


1.08E - 3 


sup 
0<t<0.1 


ip{t,x) - 


1p BC (t, X) £00 ( R ) 


9.34E - 2 


7.74E - 3 



In a second case, we alternatively choose initial data of the form 

(5.13) Vin(x) = e~ 5 ^ 2 X2 (J, sin(aO) e~ icos ^/ £ , 

i.e. 4>{x) = — cos(x). Here we find (numerically) that the caustic onset time is roughly 
given by t c « 0.24, cf. Fig. I5.il The corresponding numerical results are given in 
Fig. [KT3\ and Table\5J% 

6. Acknowledgement. The authors are grateful to Prof. Christian Ringhofer 
for fruitful discussions on this work. 
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Table 5.3 

Difference between the asymptotic solution and the Schrddinger equation for example \5.2\ for 
initial condition iPHBl (At = 10~ 4 , Ax = 1/32768J: 





e 




l 

32 


l 

1024 


sup 

0<i<0.1 


1 i>(t, x) - 




1.68E - 2 


3.19E - 3 


sup 

()<i<0.1 


ip(t,x) - 


ip sc (t,x) \\ LOO(K) 


2.73E - 1 


7.33E - 2 
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|V> cx (i, x)\ 2 , |V> cx (t, as) - ip ts (t, x)\ and |i/> cx (t, x) - ip hd (t, x)\ at t = 1.0 

Fig. 6.1. Numerical results for example \4-l\ with U(x) given by H4.2I I and & = \- We use 

-, Ax = „ * for the "exact" 



At = -rini Ax = 4n for the TS and the BD method, and At 
100 ' 64 J ' 

solution. 



100000 

A*£(t) = 5.39E - 2, Aj^(t) = 5.07E - 2, Ajf (t) = 1.56E - 2, A^ d (t) = 1.51E 



- 2. 





0.65 0.35 



(t,z) 



0.4 0.45 0.5 0.55 0.6 

ip ts (t,x)\ and |^ cx (t, x) 




Fig. 6.2. Numerical results for example \4-l\ with U(x) given by (14.21) and £ 
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Ax : 
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Ax 



for the TS, At 
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for the "exact" solution. 



Ax 



for the BD method, and At 
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131072 

A£(t) 



1.23E - 



1. A^ d (t) 



1.20E - 1, A!f (t) = 2.29E 
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|V> cx (i,x)| 2 , |V> cx (i, x) - </> ts (t,x)| and |i/> cx (i, x) - </> bd (i,x) at t = 1.0 
Fig. 6.3. Numerical results for example \4 H with U(x) given by H4.3I I and e - 



At = A , Ax = ir= for the TS and the BD method, and At = . ' , Ax = „A.„ for the 

10 62 J 1UUUUU 8192 J 

solution. 



We use 
"exact" 



A£(t) =3.47E- 



■ 3, A^ d (t) 



1.04E - 3, AJ, s (t) = 1.96E - 3, A^ d (t) = 3.65E - 4. 



ill III 12 ' ill I I 




0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.35 



|V> cx (t,x)| 2 , \i> cx (t,x) - ip ts {t,x)\ and \ip cx (t, x) - ip hd {t, x)\ at t = 0.1 
Fig. 6.4. Numerical results for example with U(x) given by 114.31 1 and e = * . . We use 



Here At 
At -- 



-, Ax 



l 



100000 

A*ci(t) 



Ax : 



for i/ie TS and At 

16384 J 

131072 f or ^ e " exac t" solution. 
, bd 



At 
100' " x 



/or t/ie BD method, and 



1.37E - 2, A^(t) = 5.52E - 3, A| s (t) = 2.76E - 



■3, A^ d (t) 



1.20E - 



22 



Z. Huang, S. Jin, P. A. Markowich, C. Sparber 




|V> ex (i,x)| 2 , \ip e *(t,x) - </> ts (t,x)| and |i/) ex (t, x) - f/> bd (*, x)\ at t = 1.0 



10 

solution. 



Fig. 6.5. Numerical results for example \4 H with U(x) given by i4Al and e 

1 A r = 



At = jn, Ax = A /or 4/ie TS and the BD method, and At ■■ 



■ . We use 
for the "exact" 



Aj£(i) = 3.26E - 2, Aj^(t) = 2.72E - 2, Af (t) = 1.51E - 2, A^ d (t) = 1.45E - 2. 



1 ...JWIIk, - 

0.3b 0.4 0.45 0.5 0.55 0.6 O.E 



1 55 0.6 0.65 




|V> ex (t,x)| 2 , \ip ex (t,x) - ip ts (t,x)\ and |^ ex (t, x) - tp (t, x)\ at t = 0.1 



Fig. 6.6. Numerical results for example \4-l\ with U(x) given by l|4.4ll and e = * . . We use 



At = 

AX : 



-, AX : 



for the TS and At = , Ax = „^ q9 for the BD method, and At 



31072 for the "exact" solution. 

A^(t) = 3.04E - 2, A^ d (t) = 4.25E - 3, A^ s (t) = 5.35E - 3, A^ d (t) = 1.21E - 
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\ip cx (t,x)\ 2 , |V> cx (i,x) - Tp ts {t,x)\ and |i/> cx (i, x) - i/j bd (t,x)\ at t = 1.0 



Fig. 6.7. Numerical results for example \4- 2\ with U(x) given by 1 14. 21 , &=\- We use At = — jjj, 
Ax = i f or the TS, At = i, Arr = ^ /or BD method, and At = ^ggggg , Ax = grog / or *' le 
"exact" solution. 

A^(t) = 3.31E - 1, Aj^(t) = 1.77E - 1, Al, s (t) = 6.16E - 2, A^ d (t) = 1.38E - 2. 




Fig. 6.8. Numerical results for example \4-2\ with U(x) given by 14.20 . e 



At = 1 

10000 

Ax - 1 



131072 
A 



Ax = 7 rS m for the TS, At 
solution. 



J_ 

id ■ 



l 

1024 ' 



Ax 



for the BD method, and At 



We use 
1 



for the "exact 



1.65, A h J(t) 



■ 9.14E - 2, ASf(t) = 2.63E - 1, A£ d (t) = 1.39E ■ 
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|?/> ex (t,x)| 2 , \ip ex (t,x) - ip ts (t,x)\ and \i/) ex (t, x) - ip bd (t, x)| at t = 1.0 
Fig. 6.9. Numerical results for example \4-^\ with U(x) given by l|4,3j l and e 



We use 



At = -^q, Ax = ^ /or f/ie TS, At = i, Ax = ^ for the BD method, and At = 1Q g QI 

3, Alf (t) = 4.02E - 2, A^ d (t) = 3.89E - 3. 



Ax = gl92 /or t/ie "exact" solution 



A%(t) = 7.30E - 2, A^ d (t) = 8.30E - 




ill 11 r I |N " ( * '^ll ljlliiii'i'i,,., 

0.4 0.45 0.5 0.55 0.1 

|V> cx (t,x)| 2 , \il> m (t,x) - il: ts (t,x)\ and |^ cx (t, x) - ip hd {t,x)\ at t = 0.1 
Fig. 6.10. Numerical results for example^. 2\ with U(x) given by 1 14. 31 1 and e = . . We use 



At = 

Ax : 



/or i/ie TS, At 



Ax 



for the "exact" solution. 



Ax 



for the BD method, and At 



A^(t) = 1.61, A^ d (t) = 9.16E - 2, Aif(t) = 2.63E - 1, A§ d (t) = 1.71E - 
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A 

V 










\ip ex (t,x)\ 2 , \ip ex (t,x) - </> ts (t,x)| and \tp c *(t,x) - i/j bd (t,x)\ at t = 1.0 

Fig. 6.11. Numerical results for example \4^\ uiith U (x) given by £0}, e = \. We use At = t^tt , 
Ax = A for the TS, At = i, Ax = A for the Bloch-decomposition method, and At — 



Ax : 



32 
1 



100000 ' 



/or the "exact" solution. 



A^(t) = 4.01E - 2, A^(t) = 5.00E - 2, A| s (t) = 1.85E - 2, A^ d (t) = 1.98E 




\ip e *(t,x)\ 2 , |i/. ex (t,x) - V ts (t,x)| and |i/> ex (t,x) 



(t, x)\ at t ; 



1 



We use 



At 
At 



Fig. 6.12. Numerical results for example \4-2\ Here U(x) is given in 114.41 1 . o — 1024 . 

Ax = 65 g 3( , /or Time-splitting, At = j^, Ax = 81 1 92 /or Bloch-decomposition, and 
• ^ x = 131072 f or ' exac t' solution. 



loooo ' 



100000 ' 

a; 



(t) = 1.35, A^(t) = 3.48E - 3, ASf (t) = 2.23E ■ 



1, A* d (t) 



1.14E - 
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J 
Ji 




0.1 


ill 


1 ; 
ill 


|V> sc Or,0.1)| 2 , 


V> sc O,0.25)| 2 


, and |^ sc (s,0.5)| 2 , e = i. 




. . Jnllllli.V . . 






L 




i 




|V> bd (:r,0.1)| : 


, |V> bd (a:,0.25)| 2 , and |V> bd 0, 0.5)| 2 , e = 




1 




1 








|</> sc O,o.i)| 2 


|</> sc (:r,0.25)| 2 , and \i> sc (x, 0.5) 


' £ 1024 ' 




1 






L 




, , , ~U 


k. 



^ bd (x,0.1)| 2 , |V> bd (z,0.25)| 2 , and \ip hd (x, 0.5)| 2 , £ = T i jf . 



Fig. 6.13. Numerical results for examvle \ 5. 2\ with U (x) given by (14. 3D . At = 10 q 00 , Arc = 32 T68 ' 
The left column shows the situation before the caustic, whereas the other two columns respectively 
present the numerical results at and after the caustic. 



